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We present a new many-parameter family of hyperbolic representations of Einstein's equations, 
which we obtain by a straightforward generalization of previously known systems. We solve the 
resulting evolution equations numerically for a Schwarzschild black hole in three spatial dimensions, 
and find that the stability of the simulation is strongly dependent on the form of the equations (i.e. 
the choice of parameters of the hyperbolic system) , independent of the numerics. For an appropriate 
range of parameters we can evolve a single 3D black hole to t ~ 600 M - 1300M, and are apparently 
limited by constraint-violating solutions of the evolution equations. We expect that our method 
should result in comparable times for evolutions of a binary black hole system. 
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I. INTRODUCTION 

A key unsolved problem in general relativity is to pro- 
vide a detailed description of the final moments of a bi- 
nary black hole system as the two black holes plunge 
together and merge. While this problem is interesting 
in its own right, the current deployment of LIGO and 
other gravitational wave interferometers provides addi- 
tional incentive for finding a timely solution: coalescing 
compact binaries are expected to be primary sources of 
gravitational radiation observable by these instruments. 
Comparison of observed gravitational wave forms to de- 
tailed theoretical predictions of binary black hole evolu- 
tion may allow one to test general relativity and other 
theories of gravitation, to identify black holes in distant 
galaxies and to measure their masses and spins. 

Although both the initial inspiral of a binary black hole 
system and the final ringdown of the resulting Kerr black 
hole are well described by perturbation theory, under- 
standing the plunge from the innermost stable quasicir- 
cular orbit through the coalescence will require numerical 
solutions of the full Einstein equations in three spatial di- 
mensions. Such numerical computations are in progress 
jl],^| ; however, they are currently plagued with instabil- 
ities that severely limit the duration of the simulations. 
Indeed, current 3D Cauchy evolution codes without built- 
in symmetries have great difficulty evolving even a sin- 
gle Schwarzschild black hole for the amount of time that 
would be required for a binary orbit. 

Many of the stability difficulties affecting black hole 
computations are undoubtedly due to the technical de- 
tails of the numerical solution scheme; there are many 
such difficulties to overcome in any large scale numerical 
solution of partial differential equations. However, there 
is also evidence that some of the stability problems are 
due to properties of the equations themselves, indepen- 
dent of any numerical approximation. In particular, by 
rewriting the equations in a different manner but leaving 



the numerical method unmodified, one can significantly 
affect the stability of the computation 

Einstein's equations, when written as a Cauchy prob- 
lem, can be decomposed into two subsystems of equa- 
tions: constraint equations that must be obeyed on each 
spacelike hypersurface, or time slice, and evolution equa- 
tions that describe how quantities propagate from one 
hypersurface to the next. An analogous decomposition 
occurs in electromagnetism, which is naturally split into 
time-independent (divergence) equations that constrain 
the fields at a particular time, and time-dependent (curl) 
equations that determine their evolution. For both elec- 
tromagnetism and gravitation, the system of equations is 
overdetermined in the following sense: if the constraint 
equations are satisfied at some initial time, then the evo- 
lution equations guarantee that they will be satisfied at 
subsequent times. For numerical black hole computa- 
tions, one typically solves the constraint equations only 
on the initial time slice, and then uses the evolution equa- 
tions to advance the solution in time. 

However, the decomposition of Einstein's equations 
into evolution equations and constraints is not unique. 
For example, one can add any combination of constraints 
to any of the evolution equations to produce a different 
decomposition. Indeed, there have been a large number 
of new formulations of 3+1 general relativity proposed 
in recent years ||,|| ^7j, many of which have attractive 
properties such as symmetric hyperbolicity. 

All such formulations must have the same physical so- 
lutions since they describe the same underlying theory. 
However, the set of evolution equations also admits un- 
physical solutions such as constraint- violating modes and 
gauge modes, and these unphysical solutions will be dif- 
ferent for each formulation. Usually one is not inter- 
ested in unphysical solutions, but if such a solution grows 
rapidly with time, any small perturbation (say, caused 
by numerical errors) that excites this solution will grow 
and eventually overwhelm the physical solution. This is 
one reason why some formulations of Einstein's equations 
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may be better suited for numerical evolution than others. 

In order to explore the extent to which different for- 
mulations of Einstein's equations affect the stability of 
numerical evolutions, we construct three new formula- 
tions of Einstein's equations, following a method similar 
to that of pi: 

1. A first-order system obtained directly from the 
ADM system. This system has five unde- 
termined constant parameters that specify con- 
straint terms to be added to the evolution equa- 
tions. These parameters determine the hyperbolic- 
ity of the evolution equations and the values of the 
characteristic speeds. We find that constraining the 
system to have physical characteristic speeds (i.e., 
the characteristic fields propagate either along the 
light cone or normal to the time slice) still leaves 
two of the five parameters free, and guarantees that 
the evolution equations are strongly hyperbolic. In 
this case, the constraint quantities also evolve in a 
strongly hyperbolic manner with physical charac- 
teristic speeds. 

2. A twelve-parameter system obtained by applying 
a parameterized change of variables to system 1. 
The additional seven parameters are completely 
free, and do not affect the hyperbolicity of either 
the evolution equations or the evolution of the con- 
straint quantities. This system can be reduced to 
either the Frittelli-Reula formulation |l3| or the 
Einstein-Christoffel formulation with an appro- 
priate choice of the parameters. The seven addi- 
tional parameters can be used either to simplify the 
equations or to improve the numerical behavior of 
the system. 

3. A two-parameter system that is obtained from sys- 
tem 2 by demanding that the principal part of the 
equations is equivalent to a scalar wave equation 
for each of the six components of <?y. This sys- 
tem is particularly simple, is symmctrizablc hyper- 
bolic with physical characteristic speeds, and in- 
cludes the Einstein-Christoffel formulation |2^] as 
a special case. 

To determine whether modifying the formulation sig- 
nificantly effects the numerical solution of the evolu- 
tion equations, we perform numerical evolutions of single 
black holes using a new 3D code we have developed. We 
evolve system 3 for simplicity. We find that by varying 
the two parameters in system 3 while keeping the numer- 
ical evolution method fixed, we can vary the run time of 
the simulation by more than an order of magnitude. For 
a single black hole, our optimum choice of parameters 
yields evolutions that run to t = 600M - 1300M. This is 
long enough that, if this result carries over to two-black- 
hole simulations, one could simulate the last few orbits 
of a binary system and the final merger. 

In section [n] we derive systems 1-3 and conditions for 
hyperbolicity. We also derive evolution equations for the 



constraint quantities and discuss their hyperbolicity. In 
section III we present numerical evolutions of system 3 for 



different choices of parameters, and show that particular 
choices yield significant improvements. In section ^ we 
discuss our results and our plans to simulate a binary 
system. 



II. PARAMETERIZED HYPERBOLIC SYSTEM 

A. 3+1 ADM 

We begin with the standard 3+1 formulation of p§| ] 
which is discussed in detail in |29f| . Four-dimensional 
spacetime is foliated by the level surfaces T, t of a function 
i(a+). Let be the unit normal vector to the hypersur- 
faces Yif Then the spacetime metric "SV„ induces the 
spatial three-metric g^ v on each S t given by 



(2.1) 



The timelike vector i M is defined such that t^t.^ = 1, 
where i ;M is the covariant derivative of t with respect to 
the spacetime metric. The lapse function N and shift 
vector (3^ are defined by 



N = -th 



(2.2) 
(2.3) 



If we adopt a coordinate system {t,x 1 } adapted to the 
spatial hypersurfaces, the line element is given in the 
usual 3+1 form 

ds 2 = -N 2 dt 2 + g ij (dx i + (3 l dt){dx 3 + (3 j dt). (2.4) 

The extrinsic curvature Kij of the spatial surfaces is given 
by 



(2.5) 



where £ denotes a Lie derivative. 

Einstein's equations are given in covariant form by 



(4 H--5 (4 V (4) ^8^, 



(2.6) 



where ^R pv and ( 4 lff are the Ricci tensor and Ricci 
scalar associated with the spacetime metric, and T^ v is 
the stress-energy tensor. In the 3+1 decomposition Ein- 
stein's equations are decomposed into the Hamiltonian 
constraint 

C = \ (R - K ab K ab + K 2 ) - 8np = 0, (2.7) 
the momentum constraints 

C, = V n X, a - V.K-SnJ, =0, (2.8) 
and the evolution equations 
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doKij = -V.-Vj-JV + NRij - 2NK ia K j a + NKK tJ 

- SnNSij - 4irN gij ( P - S), (2.9) 

where K = g ab K a b, and Vi, Rij, and R are the covari- 
ant derivative, Ricci tensor, and Ricci scalar associated 
with the spatial three- metric. The symbol do is the time 
derivative operator normal to the spatial foliation, de- 
fined by 



d = d t - £/3- 
The matter terms are defined as 



(2.10) 



Ckiij = d[kdi]ij — 0, 



(2.17) 



because second derivatives of the metric commute. 
Therefore we make the following substitution when we 
encounter second derivatives of the metric: 



dkdigij — d( k di)ij. 



(2.18) 



In terms of these new variables, the affine connection, 
Ricci tensor, and Ricci scalar are given by 



r 



klj 



-l JL I 

l (ij)k 2 u kiji 



(2.19) 



(2.11a) 
(2.11b) 
(2.11c) 



and S = g ab S a b- The definition ([D]) of the extrinsic 
curvature yields the following evolution equation for the 
spatial metric 



d ogij = 2.VA, 



(2.12) 



Note that the spatial metric and its inverse are used to 
lower and raise the indices of all spatial tensors. 



B. First-order form 



Rij — \g ah ip{id a bj) + 9ad(ij)b — 9 a dbij — d(idj^ ab ) 
I _ 1 J<1J h a rl — rl b rl a 

~r 2 tlaij 4^ U>aij ^ Ll(ij)a 2 aj bi 

+ ^d a d(ij} a + jd^ dj a b + idabj, (2.20) 



R = 9 ab 9 cd (d d d abc - d a d bcd ) + b a d a - b a b a - \d a d a 



1 rl jcab i 3 J mbc 
— 2"abc« + jO, a b c a 

The constraint equations are given by 



(2.21) 



C = \g ab g cd (d d d abc - d a d bcd ) + \b a d a - \b a b a 



— \d a d a — \d a b c d cah + ^d a b c d abc 

- \K ab K ab + \K 2 - 8tt P , 



(2.22) 



In order to cast the evolution equations in first-order 
form, we must eliminate the second derivatives of the 
spatial metric. We define a new variable (symmetric on 
its last two indices) 



dkij — dkQij , 



(2.13) 



and its traces df. = g ab dkab and bk = g ab d a bk- An evo- 
lution equation for dkij is obtained by taking a spatial 
derivative of ( 2.12| ) and using the fact that dk and <9 
commute. This yields 



dod kij = -2Nd k K tJ - 2K tJ d k N, 
where the Lie derivative of dkij is 



(2.14) 



£pdkij — f3 a d a dkij + d ai jdk/3 a + 2d ka (idj)(3 a 

+ 2g a{l d 3) d k P a . (2.15) 

Since we have introduced a new variable that we will 
evolve independently of the metric, we have an additional 
constraint, 



Ckij — dkij &k9ij 0, 



(2.16) 



which must be satisfied in order for a solution of the first- 
order evolution equations to be a solution of Einstein's 
equations. Note that the spatial derivatives of dkij must 
satisfy the constraint 



d = g ab [d a K ib - diK ab ) + \K ab d iah 
+ \K la d a -K la b a -8nJi. 



(2.23) 



Finally, the evolution equation for the extrinsic curvature 
becomes 

doKij = N Wg ah (d(id ab j) + d a d^ b — d a d b ij — 9(i4?>&) 

4- — h a rl - — rl a rl h a rl — rl b rl a 

< 2 ^aij 4 U ^aij U Uj (ij)a 2 aj bi 

+ \d a d(ij) a + jd^djab + \d a idabj — 2K ia K j a 
+KK lJ ] - didjN - ±d%d a N + d {ll) a d a N 



8-rrNSij - 4nNg tJ (p - S) 



(2.24) 



The hyperbolicity of the system of evolution equations 
can be determined by examining its principal part. Con- 
sider a system of the form 



d Q u + A l d t u = F, 



(2.25) 



where u is a column vector of the fundamental variables, 
and A 1 and F are matrices that can depend on u, but 
not on derivatives of u. One defines, for a particular unit 
one-form £j, the characteristic matrix C in the direction 
normal to 



C = A^i 



(2.26) 



The characteristic speeds in the direction & are the eigen- 
values of C. If all characteristic speeds are real, then the 
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system is said to be weakly hyperbolic. If in addition, 
C has a complete set of eigenvectors, and the matrix of 
these eigenvectors and its inverse are uniformly bounded 
functions of £j, the spacetime coordinates, and the solu- 
tion, then the system is said to be strongly hyperbolic. If 
the matrices A 1 are symmetric, the system is said to be 
symmetric hyperbolic. If the matrices A 1 can be brought 
into symmetric form by multiplying by a positive-definite 
matrix called a symmetrizer, the system is said to be sym- 
metrizable hyperbolic. Symmetric, symmetrizable, and 
strongly hyperbolic systems admit a well-posed Cauchy 
problem; weakly hyperbolic systems do not fjO[ . 

For the systems described in this paper, we explic- 
itly construct a complete set of eigenvectors that de- 
pend upon £j, the metric, and its inverse. Provided that 
the matrix norms of the metric and its inverse remain 
bounded, then the norms of the matrix of eigenvectors 
and its inverse are bounded, so the system is strongly 
hyperbolic |3l| . 

Using the method outlined in Appendix [A|, we find 
that the ADM equations written in first-order form are 
only weakly hyperbolic, as the characteristic matrix of 
the system has eigenvalues {0,±1}, but does not have a 
complete set of eigenvectors. Fortunately, the hyperbol- 
icity of the equations can be changed by "densitizing" the 
lapse and adding constraints to the evolution equations. 



C. Densitization of the lapse 



We densitize the lapse by defining 
Q = log (Ng-°) , 



(2.27) 



where g is the determinant of the three-metric, and a is 
the densitization parameter, which is an arbitrary con- 
stant. The lapse density Q and the shift vector /J 1 will 
be considered as arbitrary gauge functions independent 
of the dynamical fields. With this definition we have 



(2.28) 



diN = N {8iQ + adi) , 
didjN = N[didjQ + (diQ){djQ) + 2ad {i d j) Q 

+ ag ab d (i d j)ab - od Mh d 3 ah + a 2 d t d 3 ]- (2.29) 

Substituting the above expressions into the evolution 
equations, and examining the hyperbolicity of the mod- 
ified evolution equations, we find that densitizing the 
lapse is not sufficient to make the evolution system 
strongly hyperbolic. In order for the system to remain 
even weakly hyperbolic the densitization parameter must 
satisfy a > 0, as the eigenvalues of the characteristic ma- 
trix are now {0, ±1, ±\/2ct}. In the next section, we will 
find that densitizing the lapse is a necessary condition 
for strong hyperbolicity, and that if we demand physical 
characteristic speeds we must choose a — i . 



D. Addition of constraints: System 1 

By adding terms proportional to the constraints, we 
can modify the evolution equations for Kij and dkij with- 
out affecting the phys ical so lutio n. We modify the evo- 
lution equations (|2.14|) and (|2.24|) by 



d K l3 =(...) + lN 9lJ C + (Ng ab C a{lj)b , (2.30) 
d d ki: j =(...) + r}Ng k (iCj) + \NgijC k , (2.31) 

where (...) repre sents the right-hand side of either equa- 



tion (2.14) or (2.24), and the constraint parameters 
{l tCtVtX} are arbitrary constants. The evolution equa- 
tions are now given by 

d ogij ~ 0, (2.32) 
d Q K i:i ~ -\Ng ab [d a d bij - (1 + ()d a d {ij)b 

- (1 - C)d(id abj) + (1 + 2a)d (i d j)ab 
-19t 3 g cd d a dcdb + J9zj9 cd dad bcd ] , (2.33) 
dodkij ^ -2Nd k Kij + Ng ab (ng k{i d a K bj) + X9i 3 d a K bk 
-V9k(idj)K a b - X9z]d k K ab ) , (2.34) 

where ~ denotes equal to the principal part. For brevity, 
we have shown only the principal parts of the evolution 
equations as that is what determines the hyperbolicity of 
the system. The full evolution equations are lengthy and 
available from the authors upon request. 

We find that the eigenvalues of the characteristic ma- 
trix of the system are {0, ±1, ±ci, ±C2, ±03} where 



Cl = 


V2cr, 






C2 = 




- Ana -2% - I2ax ~ 


-3ryC, 


C3 = 




- 47 — rj — 2777 + 2\ - 


1- 47X - vC 



Thus in order for the system to be weakly hyperbolic, 
the parameters must satisfy 



cr > 0, 

n - irjcr - 2\ - I2ax - 377C > 0, 
47 - 77 - 2777 + 2% + 47x - rjC > 0. 



(2.36) 



If the above conditions are met, we find a complete set 
of eigenvectors, so that the system is strongly hyperbolic, 
unless one of the following conditions occur: 



Cj = 0, 

ci = c 3 ^ 1, 

Cl = c 3 = 1 7^ c 2 . 



(2.37a) 
(2.37b) 
(2.37c) 



If any of the above conditions are met, the system is 
only weakly hyperbolic. Note that if a = 0, then ci = 0, 
so that densitizing the lapse is a necessary condition for 
strong hyperbolicty. Also note that if 77 = \ = 0> then 
c 2 = 0, so that constraints must be added to the evolution 
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equation for d k ij in order to have a strongly hyperbolic 
system. 

For physical characteristic speeds, each of the c, is ei- 
ther zero or unity. To make them all unity (the only 
choice that yields strongly hyperbolic evolution equa- 
tions) requires cither 



a = 1/2, 

8 + 57/ + 107?y 



c 



77(7 + 67) 
4 + 67 — 77 — 3777 
(7 + 67) ' 



(2.38a) 
(2.38b) 

(2.38c) 



{<j,7,C,f/-,x} = {i,-|,-|(23 + 20x),|,x}- 



In the first case, there are two free parameters, and in 
the second case there is one. In both cases, the evolu- 
tion equations are strongly hyperbolic as long as the free 
parameters are chosen such that all five parameters are 
finite. 



E. Evolution of the constraints 

Taking do of the constraints, and replacing all deriva- 
tives of the fundamental variables with the constraints 
and their spatial derivatives, we obtain the following 
equations for the evolution of the constraints: 

doC^-^(2- V + 2 X )NgPW p C q , (2.40) 
dod ~ -(1 + 2 1 )Nd l C + \Ng^g rs [(1 - Qd q C prsl 

+ (1 + Qd p C slqr - (1 + 2a)d p C qirs ] , (2.41) 
d C kij ~ 0, (2.42) 

doCkUj = \r]N (gj[id k ]Ci + 9i[id k ]Cj) 

+ X Ng ij d [k C l] , (2.43) 

where again for brevity we have only shown the principal 
parts of the equations. 

The eigenvalues for the constraint evolution system are 
{0, ±c 2 , ±c 3 }. Because this is a subset of the eigenvalues 
of the evolution equations, the constraints will propagate 
at the same speeds as some of the characteristic fields of 
the evolved quantities. Furthermore, we find that the 
constraint evolution system is strongly hyperbolic when- 
ever the regular evolution system is strongly hyperbolic. 

F. Redefining the variables: System 2 

The evolution equations can also be modified by re- 
defining the variables that are evolved. We define the 
generalized extrinsic curvature Pij using the relation 



where z is an arbitrary parameter. The inverse transfor- 
mation is given by 



where P = g P a b, and 



1 + 31' 



(2.45) 



(2.46) 



which implies that z ^ — ^ for the inverse transformation 
to exist. 

We define the generalized derivative of the metric, 
M ki j, using the relation 



(2.39) M ^3 = § [kdkij + ed (lf)k + g l3 

+g k (i cdj) + dbj) I . 



ad k + bbi. 



(2.47) 



The inverse transformation is given by 



Pi, 



Ki. 



zg lj K, 



(2.44) 



d Hj = 2 {kM Hj + eM {m + g l3 [dM k + bW k ] 

+g k{l [cM 3) + dW f) ] } , (2.48) 

where the traces M k = g ab M ka i, and W k = g ab M a t, k , and 

5a = 6bce — 6dde — ae 2 + be 2 + ci 2 — de 2 + 8bck 

- 8ddk - Adek + 2bek + 2cek - Aak 2 , (2.49a) 
6b = -8bce + 8dde + 2de 2 - 2ce 2 - 4bck + 4adk 

+ Aaek - 2bek + 2dek - 4bk 2 , (2.49b) 
6c = -8bce + 8ade + 2ae 2 - 2be 2 - Abck + 4adk 

+ Adek - 2cek + 2dek - Ack 2 , (2.49c) 
6d = Abce - 4ade - Aae 2 + \2bck - \2adk + Abek 

+ Acek~Adk 2 , (2.496) 
6 a e = 2e, (2.49e) 
S k = -i-2k, (2.49f) 
6 = i 2 - ek- 2k 2 , (2.49g) 
5 = <5 (10fcc - lOad -ae + 3be + 3ce + de + e 2 

- 6dk - 2bk - 2ck - Adk -ek- 2k 2 ). (2.49h) 

For the inverse transformation to exist, 5 7^ 0. 

Thus we have seven additional redefinition parameters 
{a, b, c, d, e, fc, £} (or equivalently {a,b,c,d,e,k,z}) that 
can be used to mod ify t he ev olution equations. Note 
that equations (2.46) and ( 2.49| ) remain true under inter- 
change of {a, 6, c, d, e, k, z} and {a, b, c, d, e, k, z}. 

When the principal terms in system 1 are transformed, 
terms containing deriv atives of th e metric appear because 
of the traces in (2.45 ) and ( 2.48 ). These terms are elim- 
inated using (2.12) and ( [2.15 ). 

The redefinition parameters do not change the eigen- 
values of the evolution system, nor do they change 
whether or not the system is strongly hyperbolic (see 
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Appendix ||). In addition, they have no effect on the 
principal part of the constraint evolution equations. The 
redefinition parameters, however, do affect the eigenvec- 
tors of the evolution system and thus they also affect the 
characteristic fields. In addition, the redefinition parame- 
ters change the nonlinear terms in the nonprincipal parts 
of the evolution equations and the constraint evolution 
system. 

The principal parts of the evolution equations for 
and Mkij are 

dogij 0, (2.50) 
d Pi 3 ~ -Ng ab (fiidaMbij + iJ, 2 d a M^ b + ^s,d^M abj) 
+ fiid^Mj) ab + ^ 5 gijg cd d a M cdb 
+^g lJ g cd d a M bcd ) , (2.51) 
d M k ij ^ -N {vid k P i:j + v 2 d {i P j)k + v 3 g al 'g k{i d a P bj) 
+ v±g l3 g ah d a P bk + ^55 ab 5fe( i a j )Pah 
+is 6 g t3 g ab d k P ab ) , (2.52) 

where 

Hx = k - |(1 + C)e, (2.53a) 
M 2 = i(l-C)e-(l + 0^ (2.53b) 
Ms = (1 + 6cr)6 - (1 - ()k - |(1 - 4cr - 3C)d 

+ |(1 +4o- + C)e, (2.53c) 
M4 = (1 + 6cr)a + (1 + 2a)k - i(l - 4a- - 3()c 

-i(l-C)e, (2.53d) 
/j 5 = (1 + 27 + Az + 672 + 6az)b - (7 + 2z + 3jz)k 

- |(1 + 27 + Az + 67Z - Aaz + C)c 

+ i( 7 + 2z + 37z + 4cjz)e, (2.53e) 
/i6 = (1 + 2 7 + 4z + 67Z + 6az)a 
+ (7 + 2z + 372 + 2az)k 

- |(1 + 27 + 4z + 67Z - 4(tz + C)d 

- i( 7 + 2z + 37i)e, (2.53f) 

vi = k, (2.53g) 

v 2 = e, (2.53h) 
1/3 = |(2 - 277 - X )d- §(77 + 3 X )c - ±(77 + 2 X )e 

- \r]k, (2.531) 
^4 = |(2 - 27? - X )6 - |(77 + 3 X )a - ±7?e - i*fc, (2.53j) 
^5 = I (2 + 77 + 3x + 6z + 2r)z + 6x^)c 

+ i(2?7 + x + 2z + 4?7z + 2xz)d + ±(77 + 2rjz)k 

+ ±(r7-t-2x + 4z + 2r7z + 4x^)e, (2.53k) 

1/6 = |(2 + 77 + 3x + 6z + 27yz + 6xz)a 

+ |(2jj + x + 2z + 4?7z + 2xz)6 

+ ±(77 + 2772)13 + i(x + 2z + 2xz)fc. (2.531) 

Again, the full evolution equations are available from the 
authors upon request. 



Furthermore, we note that if (i{ = nvi for all i and 
constant k, the system is symmetrizable hyperbolic us- 
ing the energy norm argument of JT3J . These conditions, 
however, do not have to be met for the system to be well- 
posed. It is possible to construct a symmetrizer for any 
of the strongly hyperbolic systems. 



G. Evolving with contravariant indices 

So far, we have written all of our fundamental variables 
with covariant indices. Alternatively, we could have de- 
fined the new variable 



dug 1 - 



(2.54) 



Note that d k i 



-D 



kij- 



If we evolve {g ij , P lj , M k lJ } 



instead of {gij, Pij , Mkij}, it would result in only triv- 
ial changes to the principal parts of the equations. The 
characteristic speeds would be unchanged, as would the 
nature of the hyperbolicity of the system, since the prin- 
cipal part of the metric evolution equation is zero (See 
Appendix [^) . The only changes would occur in the non- 
linear terms of the evolution equations. 



H. Frittelli-Reula system 



We recover the system of 13 1 if we make the following 
choices for our parameters: 



W,~fX,v,x} = {- 



;l+3a 27 



- 1 4 

1+3/3' ' ' 1+3 



{z, k, a, b, c, d, e} — {f3, 1, a, 0, 0, 0, 0}, 



%}, (2.55a) 
(2.55b) 



where {a, (3, 7, e} correspond to {a, f3, 7, e} in [ p^[ . How- 
ever, as pointed out in [|T|, this system is not symmetric 



hyperbolic unless the term —2h l ^M 3 ' z k in Eq. (16) of 
fbl is replaced with —2h L ( l M J ' k k , by adding a term pro- 



portional to the constraint (2.17). In our syste m this 
corresponds to changing £ = 1 to £ = — lin ( 2.55[) . 

In P j3l| , p2| |, this correction has been made for the pa- 
rameter choice {a,/3, 7 , e} = { — 1,-1,1,^}; we recover 
this system if we choose our parameters to be 

{a, 7 ,C,»7,x} = {i-l,-l,4,-2}, (2.56a) 
{z, k, d, b, c, d, e} = {-1, 1,-1,0,0, 0, 0}. (2.56b) 



The system of 31 3^] was further generalized by |^7) , 
who used the constraints to modify th e evo lution equa- 
tions in a manner similar to that in Sec. [ID . We recover 
the system of |27| by choosing: 

W, 7, C, V, x} = {i -7, 26 - 1, 4r7, -2r)}, (2.57a) 
{z, k, a, b, c, d, e} = {-1, 1, -1, 0, 0, 0, 0}, (2.57b) 

where {7, 0, 77} correspond to {7, 6, 77} in pTj. 
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I. Einstein-Christoffel system 



III. NUMERICAL RESULTS 



We recover the system of |22| if we make the following 
choices for our parameters: 

7>C.»»,X} = {§,0,-1,4,0}, (2.58a) 
{z, k, a, b, c, d, e} = {0, 1, 0, 0, 2, -2, 0}. (2.58b) 

This system is symmctrizable hyperbolic and has a very 
simple principal part 



doPij 
d M ki j 



-Ng ab d a M bij , 
-Nd k P l3 . 



(2.59a) 
(2.59b) 



Essentially this system is a set of six (one for each {i,j} 
pair) coupled quasilinear scalar wave equations with non- 
linear source terms. 



J. Generalized Einstein-Christoffel: System 3 

If we examine the principal part of System 2, and de- 
mand that \x\ = v\ = \ and all other /x, and Vi vanish, we 
obtain a two-parameter system {77, z} that has the same 
simple wave-like form (2.59) as the Einstein-Christoffel 
system. This system is obtained by setting 



{a, 7, C, rj, X } = a, -1, V, (2.60a) 



{z, k, a, b, c, d, e} = {z, 1, 4+v 2 1 ^ z+9r ' z 



2 ' 217 
-4 

1 3±"j=ia! 2,-2,0}, 



2,i 



(2.60b) 



where z^-^ and rj 7^ 0. This system has physical char- 
acteristic speeds and is symmctrizable hyperbolic. The 
free parameter 77 will affect the principal part of the con- 
straint evolution equations, while the parameter z will 
affect only the nonlinear terms in the evolution equa- 
tions and the constraint evolution equations. It is this 
system that we will explore numerically in the follow- 
ing section. The complete equations for this system are 
available upon request from the authors. 

The characteristic eigenfields of this syst em a re partic- 
ularly simple, and can be obtained from (2.59) without 
the use of the lengthy decomposition procedure described 
in Appendix^. In a direction £j, the eigenfields are 



U° 



± _ 



p ij ±eM kiJ . 



(2.61a) 
(2.61b) 
(2.61c) 



The U° quantities propagate along the normal to the 
time slice (coordinate speed — f3 l ), and the quantities 
propagate along the light cone (coordinate speed — (3 l ± 
NC). 



In this section we present results from a numerical code 
that solves the evolution equations of System 3 in three 
spatial dimensions plus time. This code, which will be 
described in detail elsewhere [[33| , is a three-dimensional 
generalization of a spherically symmetric code discussed 
previously H|, and is based on pseudospectral colloca- 
tion methods. Our code works in full 3D; we do not 
exploit any symmetries of the black hole solutions that 
we evolve. 

In this paper, we will concern ourselves only with single 
black hole spacetimes. In this case, we solve the evolu- 
tion equations in a spherical shell extending from inside 
the horizon to some artificial outer boundary. Although 
we use standard spherical polar coordinates (r,9,<ft), we 
evolve the Cartesian components of our variables; this 
allows us to use scalar spherical harmonics Yf. m {6, <f>) as 
angular basis functions for all quantities. We use Cheby- 
shev polynomials as the basis functions in radius. 

As described in jj4|, we use the method of lines in 
order to integrate forward in time with a fourth-order 
Runge-Kutta method. Boundary conditions are imposed 
by constructing the characteristic fields that propagate 
normal to the boundary, and imposing conditions only 
on those fields that propagate into the computational do- 
main. Since all characteristic fields at the inner boundary 
are outgoing (into the hole), no boundary condition is 
needed there and none is imposed. At the outer bound- 
ary, we impose dtU~ = on each of the characteristic 
fields U~ that is ingoing there. We use analytic initial 
data corresponding to time-independent slicings of a sin- 
gle black hole, and fix the gauge quantities Q and (3 l to 
their analytic values. Note that the constraint equations 
are not solved explicitly, but are instead used as a check 
on the accuracy of our numerical integrations. 



A. Einstein-Christoffel 

Figure |l| shows the £2 norm of a component of 
the momentum constraint for several evolutions of a 
Schwarzschild black hole using the Einstein-Christoffel 
system, which is equivalent to System 3 with rj = 4 and 
z = 0. Initially the fields are given analytically on a 
Painleve-Gullstrand time slice [p5|-p8| . Explicit formulae 
for our variables on the initial slice can be found in p4| . 

As is evident from Figure [j], the constraint increases 
with time until the simulation terminates. The evolu- 
tions with higher radial resolution run longer, but in- 
crease at approximately the same rate. In addition, for a 
fixed resolution, we see no significant dependence on At, 
and for a fixed radial resolution and time step, we see no 
significant dependence on angular resolution. This sug- 
gests that the growth of the constraints may be due to an 
unphysical solution of the equations rather than a numer- 
ical instability. Numerical instabilities typically become 
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FIG. 1. Momentum constraint C x versus time for evolu- 
tions of a Painleve-Gullstrand time slicing of a Schwarzschild 
black hole using the Einstein-Christoffel system. Results are 
plotted for several radial resolutions ranging from N r — 10 to 
N r — 40, fixed angular resolution £ — 7 and fixed time reso- 
lution At = 0.015M. Higher radial resolutions correspond to 
smaller errors. 

worse when one increases the resolution or decreases the 
time step. In contrast, our results appear consistent with 
an unphysical solution of the equations that initially has 
a nonzero amplitude because of small numerical errors. 

B. Generalized Einstein-Christoffel 

Because we suspected that the instability shown in Fig- 
ure [l] is related to the equations rather than the numer- 
ical method, we repeated the above evolutions for vari- 
ous values of the free parameters r\ and z, searching the 
two-dimensional parameter space for systems of evolution 
equations that might be better-behaved. We found that 
for 77 ~ 4/33 and z ~ —1/4, our numerical simulations 
ran for an order of magnitude longer than for the basic 
Einstein-Christoffel system. Typical results are plotted 
in Figure^. Although a growing mode is still present, its 
growth rate is much smaller than in Figure [ij and the mo- 
mentum constraint is less than 10~ 3 until approximately 
600M. 

We see no evidence that the growth is due to a numer- 
ical instability. In contrast, the evolutions in Figure ^ 
appear to converge to a well-defined solution. This so- 
lution is the sum of two components: a roughly time- 
independent component and an exponentially growing 
component. By extrapolating backwards along the grow- 




200 400 600 

t/M 

FIG. 2. Momentum constraint C x versus time for the same 
evolutions shown in Figure |l| except r\ = 4/33 and z = —1/4, 
and we plot more radial resolutions. If the outer boundary is 
moved out to r = 40M, the run time extends to ~ 1300M for 
the same accuracy. 

ing component in Figure g one can see that this compo- 
nent has mag nitude ~ 10" 16 at t = 0, which is on the 
order of machine roundoff error. 

As in the Einstein-Christoffel case, we see no depen- 
dence on angular resolution or on At. Our results do de- 
pend upon the location of the outer boundary. In the evo- 
lution shown in Figure ||, the spherical domain extends 
from r = 1.9Af to r = 11. 9M. Moving the outer bound- 
ary further out results in longer evolutions, increasing 
the run time from around 600M up to 1300M with the 
outer boundary at r = 40M. Moving the outer bound- 
ary beyond r — 40M, however, does not have any effect. 
We suspect that errors in the constraints at the outer 
boundary may be responsible for the constraint-violating 
modes. 

In addition to Painleve-Gullstrand slicings, we have 
run Kerr-Schild [ ^9|j40| and harmonic-time jyJflD slic- 
ings of a Schwarzschild black hole with similar qualitative 
results. For example, using the parameters of Figure ^ 
with a Kerr-Schild slicing as initial data, we were able 
to evolve up to t = 500M with the outer boundary at 
r = 11.9M, and up to t = 90(W with the outer bound- 
ary at r = 40M. We have also evolved a Kerr black 
hole with a = M /2 to t = 400M, with a spherical shell 
extending from r = 1.5M to r = 11. 5M. 



IV. DISCUSSION 

We have constructed a twelve-parameter family of 
hyperbolic formulations of Einstein's equations that is 
strongly hyperbolic for a wide range of the parameter 
space, and that includes the systems of Jl3|] and p2fl . By 
restricting ourselves to a two-parameter subset of these 
equations, we have demonstrated how the choice of pa- 
rameters can have a dramatic effect upon the amount 
of time a numerical simulation of a black hole can run 
before being swamped by an unphysical solution. 

Our runs with our best parameter choices appear to be 
limited only by the growth of constraint- violating modes 
which grow from the level of numerical roundoff errors. 
At present, we have no explanation as to why the par- 
ticular choice of parameters used to produce Figure || 
is so much better than the Einstein-Christoffel system. 
This choice was found empirically by running our code 
for various values of the parameters. It would be ex- 
tremely useful to have some theoretical understanding of 
why one particular parameter choice behaves much bet- 
ter than another, as the cost of performing a parameter 
search on the full twelve-parameter system would be pro- 
hibitive. 

Having found a system of equations and a numerical 
method capable of evolving a single black hole for a phys- 
ically interesting length of time, we now plan to turn 
our attention to the evolution of a binary black hole sys- 
tem. We expect our computational method to be capable 
of evolving the binary system to times on the order of 
several hundred M given appropriate gauge conditions. 
When we realize this, we will be able to simulate the last 
orbit or two prior to the plunge as well as the coalescence 
itself. 



solve C'wi = XiWi- The eigenvalues of the original matrix 
C are A^, and the eigenvectors are D~ 1 u>i. 

The transformation D is the decomposition of each of 
the fundamental tensor (or tensor-like) quantities into its 
irreducible parts, as we now describe. Suppose v = Du. 
Then if u and v are scalars, D is the identity operator, 
v = u. For a vector quantity u — Vi, D is defined by 



V i = D- 1 v = vf r) +&V( L \ 



(Al) 



where the longitudinal and transverse parts of V. are 
given by 



y(T) _ , to y 
i L i v m, 



(A2a) 
(A2b) 

where _Ly is the projection operator 

Uj=9ij-^3- ( A3 ) 
For a symmetric second-rank tensor u = Pij , 



Pij = D^v 



_ p(TTs) 

ij 



+ k(9ij-tei)P, 



(A4) 



where 



p 


— mn p 

— y r mn i 


(A5a) 


p{LL) 


— {men p 

— S S r mn j 


(A5b) 


i 


_ trn 1 n p 

— S r mnt 


(A5c) 


p(TTs) 

ij 


/ i mi n 1 | | mn\ p 

= ^ (« J') _ 2 - L y'- L J r mn- 


(A5d) 
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APPENDIX A: HYPERBOLICITY 

To determine the characteri stic speeds and eigenvec- 
tors of a system of the form Q2.25 ), we proceed in two 
steps. Instead of directly finding the eigenvalues and 
eigenvectors of C = A 1 ^, we first construct a transfor- 
mation D such that C = DA l £iD~ 1 is independent of 
the direction and of the metric quantities Qij . We then 



For a third-rank object u — M k ij, symmetric on its 
last two indices, 



M kij =D~ 1 v 



= M, 



(TTT) 
kij 



+ \m£ ll) ( 9j)k -^ k ) 
+ \mI llt) Onj-^j) 

+ M[f LT) (3^ k -g j)k ) 

+ i M ™ fe)a _ 5 . )fe ) 
+ i M (w) & ._ gij) 

+ 1 M (LRR) fa gij - fa^) 

+ M( RR V(g kii Z J) -Z k tiZ J ), (A6) 
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where 



M (RRL) 


— 9 S lvl cabi 


(A7a) 


M (LRR) 


= q ab £ c M cab , 


(A7b) 


M (LLL) 


= ec^ b M cab , 


(A7c) 


-Jyr(RRT) 

lvl i 


= g ca ± i b M cabj 


(A7d) 


j^(TRR) 


= g ab ± % c M cab , 


(A7e) 


M (LLT) 




(A7f) 


M (TLL) 


= re 6 v M cab , 


(A7g) 




= e (± l a ± b -\ u^ ab ) M cab , 


(A7h) 


M {TTLs) 


= 4 b (vv -| u 3 ±- ca ) M cab , 


(A7i) 


M (TTLa) 

^3 


= £ - i -u c - i -j] a M cab , 

L*" J I 


(A7j) 


kij 


= \± h c ± 4 a ±} 

I k i j 






1 1 /q 1 ab 1 c 1 ca i f) 






- | J_ fcj (3 T cb T 4 a - J."^ 6 - 


± ab ± l c ) 




- \ ± ki (3 1™!./ - ± ab ±f 






~ ±° b ±f)] M cab . 


(A7k) 



Finally for a four-index object , symmetric on its 
last two indices and antisymmetric on its first two indices, 



Ckiij =D v 



r (TF) 3 r (TTRRa) 4 r (RTTRa) 
U klij 5 U kl 9ij + 5 U kl 

. 2 fMTTRRa) n (TTRRa) 



9ij 



5 \ i[k 

6 / q{RTTRo) ^ 



. n {RTTRa) 
5\^i[k 9l]3+ L j[k 91] 

-icir%9^icir R) 9^ 



i r {RLTR) 



S^q(RLTR) 



, 28 n (RLTR) e 

9l](m + Ts C (i 9j)[kK,l] 



^ 15°[fc 

+ 5 L [k &]9ij ~ TE L [k 9i](M 

^ TLR) 93)[k ^ + 2^ 9m ^C( R ^\ (A8) 



where 



q(RLLR) 


— n cb(d(ap 
= 9 t ? L-cdab, 


(A9a) 


MRTLR) 


— „cb£a | d n 

= 9 4 -Lj <-<cdab, 


(A9b) 


p(RLTR) 


— „cb£d | a 

= .9 4 -L» <--cdab, 


(A9c) 


p(LTRR) 


= 5 4 - L l <-- C(iab , 


(A9d) 


p(RTT Rs) 
U ij 


— ^cfc / I I a 1 | \ da\ n 

= 9 -J-j-) ~2 -Ly"- 1 ) L cdab, 


(A9e) 


r (RTTRa) 


— n cb 1 d 1 a ^> 
= 5 - L [i -Lj] <-cciafc, 


(A9f) 




— „afc | ci d ^ 

= 5 - L [i - L j] ^-odabj 


(A9g) 



t(TF) 
'klij 



9k C 9i d 9i a 9 j b 



28 „cba d 
JE9 9 (i9j)[k9l] 



^9 cb 9 d {i9j)[k9i ] 
■l9v9 ab 9k9i d )C c dab. 



ffl°V (j 5i)[fcSq' 



(A9h) 



Strictly speaking, eqs. ( |A6|) and ( |A8[) are not complete 
irreducible decompositions. However, they are sufficient 
for our purposes. 

If u consists of several tensor (or tensor-like) objects, 
then the effect of D is to transform each object inde- 
pendently according to the above definitions. In matrix 
language, this means that D is block diagonal. 



APPENDIX B: CHANGE OF VARIABLES AND 
HYPERBOLICITY 

In this section we show that for a system of the 
form (2.25), a change of variables (such as the trans- 



formation from system 1 to system 2, or the raising and 
lowering of tensor indices of fundamental variables) does 
not change either the characteristic speeds or whether 
the system is strongly hyperbolic, provided that the fol- 
lowing conditions are met: 

1. The change of variables is linear in all dynamical 
variables except possibly the metric. 

2. The change of variables is invertible. 

3. Time and space derivatives of the metric can be 
written as a sum of only non-principal terms (for 
example, using (2.12) and (2.13)). 



For a system of the form ( 2.25 ), we choose an arbi- 
trar y dir ection £i and we define the matrix C according 
to ( 2.26 ). The system has k characteristic speeds \( k > 
and eigenvectors uA fe ) that obey 



(Bl) 



If M is the matrix whose columns are the eigenvectors 
u/ fe ', then strong hyperbolicity is equivalent to detM ^ 
0, with all real. 

Now consider a change of varia bles v — Tu, where T 
is a matrix. If we multiply ( |2.25 ) on the left by T, we 
obtain 

d v + TA l T- 1 d l v = TF+ (d T)u + TA l T- 1 (d l T)u 

= F' . (B2) 

In the last step, we have used property 1 above to rewrite 
d{F and doT in terms of derivatives of the metric, and 
we have used property 3 to eliminate these derivatives, 
absorbing the resulting non-principal terms into the new 
right-hand side F 1 . 



The characteristic matrix for (B2) in the direction 
is C = TA l T~ 1 £ y . l . Note that 
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C'Tw (k) = TA'T^^Tw^ = TA l £, iW {k) = X^Tw (k \ 

(B3) 



so @ and (|2725| ) have the sam e characteristic speeds 
\( k \ and the eigenvectors of (B2) are Tw^ k \ 



Furthermore, the matrix of eigenvectors for (B2) is 
M' = (TM) T , so 



detM' = det(TM) T = det T detM. 



(B4) 



If the transformation T is invertible, det M' 7^ if and 
only if det M 7^ , so (B2) is strongly hyperbolic if and 
only if (2.25) is. 
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